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Abstract . Recent work on analytical satellite- 
perturbation theory has involved the completion of 
a revision to 4th order for zonal harmonics, the 
addition of a treatment for ocean tides, an exten- 
sion of the treatment for the noninertial reference 
system, and the completion of a theory for direct 
solar- radiation pressure and earth-albedo pressure. 
Combined with a theory for tesseral-harmonics , 
lunisolar, and body- tide perturbations, these 
formulations provide a comprehensive orbit- computa- 
tion program. Detailed comparisons with numerical 
integration and observations are presented to as- 
sess the accuracy of each theoretical development. 

Introduction 

Attempts to find analytical descriptions of 
satellite motion predate the age of artificial 
earth satellites. The considerable body of theory 
in existence at that time became the foundation on 
which to build solutions to specific problems 
arising from the desire to calculate trajectories 
of artificial satellites. The celebrated volume 
64 of the Astronomical Journal can be considered 
the beginning of the field of celestial mechanics 
for satellite geodesy. Three papers in particular 
appeared in that issue [Brouwer, 1959; Garfinkel, 
1959; Kozai, 1959]; of these authors, Garfinkel 
and Kozai are contributing to the field today. 

Those articles all address essentially the same 
problem — perturbations due to Jj, J 3 , and J 4 — 
and are significant in two respects; That problem 
is still receiving attention, and the methods em- 
ployed then are still in use. One part of this 
paper is devoted to what is called the main problem 
of satellite theory; some of the present results 
will be reviewed. The method used then by Kozai 
to integrate the Lagrange planetary equations is 
almost commonplace now. The device of canonical 
transformations employed by Brouwer and Garfinkel 
is now used almost exclusively when higher order 
solutions are developed for specific problems. The 
von Ziepei method of finding a canonical trans- 
formation led to a significant generalization by 
Hori [1966, 1973] in the method of Lie Series, a 
method that has become the sine qua non of modern 
methods . Just because the beginning work was 
similar to the present, however, does not mean 
that no progress has been made . On the contrary , 
enormous strides have been taken and considerable 
work probably remains . 

The main motivation for developing elaborate 
analytical descriptions for satellite motion is to 
aid in understanding the forces causing the motion. 
A second practical reason is the potential economy 
available for certain applications. As our under- 
standing of the driving forces increases and as 
the observational accuracy improves, the require- 
ments for accuracy become correspondingly strin- 
gent. An accuracy goal — say 1 cm — is easily 
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set, but how to verify that an accuracy has been 
achieved is not so clear. 

In 1967, I presented a paper similar to this 
one [Gaposchkin, 1968] , describing a philosophy 
of how to develop, combine, and verify a complete 
satellite theory and outlining the status at that 
time. Basically, the method consisted of two 
steps — a comparison with numerical integration to 
verify that the mathematical problem had been 
properly solved, followed by a comparison with ob- 
servations to verify that the mathematical problem 
was an adequate description of the physical prob- 
lem. I am not so optimistic now. First, preserv- 
ing the accuracy of numerical-integration methods 
for long time periods poses considerable problems 
[Velez, 1975; Balmino, 1975] . In this context, a 
long period is measured in terms of the number of 
revolutions of the fastest body, normally the 
satellite; integrations for more than 1000 revolu- 
tions of anything are difficult and time consuming. 
Although the theory can be used to test the inte- 
gration, rather than the other way around, even 
the use of numerical integration to test short- 
period perturbations has proved difficult. Second, 
verifying the theory by means of data analysis 
presents significant problems. Many cases occur 
in which several forces have similar qualitative 
orbital effects that tend to cancel. A good 
example is shown by lunar and solar perturbations, 
where the direct effects, the effects due to tidal 
deformation, and the effects due to precession and 
nutation all have the same origin and produce 
perturbations with the same spectral character; 
in some cases, these effects add, and in others, 
they cancel. A third aspect of this verification 
process is the necessity to know certain physical 
quantities. While this is not a problem in the 
comparison with numerical integration, it is 
critical in the analysis of observations . In 
reality, then, the verification must be combined 
with the determination of physical quantities. 
Finally, although individual components of the 
satellite theory can be verified, in practical 
terms, the theory must all fit together and sev- 
eral interactions should be taken into account. 
Indeed, there are significant difficulties, but 
the situation is far from hopeless, and in the 
following, I describe where we are today in the 
theory, the verification, and the data analysis. 

Selection of Variables 

The first thing to select is a set of variables 
to be used for the analysis. The most popular set 
is the Kepler elements (ou = argument of perigee, 

S 2 = right ascension of the ascending node, I = in- 
clination, e = eccentricity, M = mean anomaly, and 
a = semimajor axis) . In practice, a is obtained 
from the mean motion n according to Kepler's third 
law n a J = constant, and n then becomes the sixth 
variable. This is done for the practical reason 
that n is the more easily and more accurately 
determined quantity, and a becomes a derived 
quantity. This set of variables has the conceptual 
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advantage that each has a simple physical meaning. 
Kepler elements have the one drawback that for 
zero eccentricity or inclination, they become 
degenerate. This degeneracy is also a problem for 
small eccentricities and inclinations, in that 
these variables become highly correlated in any 
adjustment procedure that attempts to determine 
them by a statistical process using observations. 
Finally, this degeneracy presents analytical dif- 
ficulties; some series expansions become much 
longer than they would if another set of variables 
were chosen. Nevertheless, Kepler elements con- 
tinue to be the most widely used both for analyti- 
cal work and for reporting results . 

For the reason cited above, some fundamental 
analysis is now done with another set of variables; 
a list of candidate variables was given by Gaposch- 
kin [1973] . For the following, the inclination 
degeneracy has not been addressed, but the more 
important and relevant degeneracy in eccentricity 
has been overcome by using the nonsingular vari- 
ables 

£ = e cos u, n = e sin co, M + u, fi, 

2 

L = a, H = G cos X , 

where 

G 2 = L 2 (l - e 2 ) = L 2 (l - £ 2 - n 2 ) . 

In that set of variables,* the Kepler elements are 
recognizable, along with the Delaunay variables 
l = M, g = ui, h = ft, L, G, and H. The Delaunay 
variables are used to derive the long-period 
and secular perturbations and are then formally 
combined into the set of nonsingular variables 
for numerical evaluation. It can be shown that 
these nonsingular variables satisfy the d'Alembert 
characteristic with respect to eccentricity, al- 
though not with respect to sin I. The d'Alembert 
characteristic is that the .lowest powers of e and 
sin I of the coefficient (kM + qw) in the 

trigonometric series are a and 3, respectively, 
where 

a = |k - q| (mod. 2) 

3 = | q | (mod. 2) 

Therefore, any expression satisfying the 
d'Alembert characteristic is well behaved as e 
(or sin I) goes to zero. If it is necessary to 
compute perturbations for small inclinations , then 
the variable M + w + satisfies the d'Alembert 
characteristic and is suitable. 

For the main problem, we have formally obtained 
expressions for the nonsingular variables analyt- 
ically. However, if the computer word has 
sufficient accuracy, then the nonsingular variables 
can be calculated numerically, and the well- 
conditioned properties of these variables can be 
realized. Therefore, it is entirely adequate, 
if necessary, to derive perturbations in Kepler 
elements and numerically combine then into non- 


* Expressions are often written in this expanded 
set of variables, but formally, this is only a 
notational convenience. 


singular variables for calculating an ephemeris . 

Of course, a unified treatment in nonsingular 
variables would be preferable, although it is not 
always the most convenient solution. We can then 
develop perturbations in the most convenient set 
of variables for the particular problem and then 
unify the variables at the calculation stage . 

Mean elements are a key to the construction of 
an analytical theory. In the framework of per- 
turbation theory, the mean elements are the zero- 
order reference for the development. They become 
the constants of integration and therefore play 
a similar role to the initial conditions in solv- 
ing differential equations . Each perturbation 
theory has implicit in it a definition of mean 
elements, and generally the relation between mean 
elements and the initial conditions does not 
receive any attention. 

In the present situation, several perturbation 
theories are employed in the same computation , 
with the mean elements being empirically obtained 
from observations. This situation is rigorously 
correct when a single perturbation theory is 
employed, provided sui table partial derivatives 
are available. In the general case, the mean 
elements must have the same formal definition, to 
the accuracy of the theory. 

A second aspect of mean elements concerns their 
constancy. If they were truly constants of the 
motion, if we knew all the numerical constants 
entering the theory, and if our observations were 
without error, then the mean elements for a 
satellite would be the same at different epochs. 

Any variations in them would have to be due to 
errors in the theory, errors in some numerical 
constants, or errors in the data. Assuming that 
we can control errors in the theory and the data, 
then the variations in the mean elements can be 
used to get information on the numerical constants 
(i.e., the physical parameters) entering the orbit 
theory. In fact, this has been the basis of much 
of the geodetic and geophysical information 
obtained from satellite data. 

Analytical Methods 

Basically, three methods are used to develop 
satellite theory. To begin with, we cannot hope 
to find exact closed-form solutions to the equa- 
tions of motion and must seek approximate solutions 
by some perturbation method. The simplest method 
is to recast the equations of motion in our chosen 
set of variables, for example , Kepler elements. 

This results in a rigorously equivalent set of 
six coupled first-order differential equations, 
called the Lagrange planetary equations. For 
some perturbations (for example, for tesseral 
harmonics), we can expand these variables around 
a reference orbit, say a precessing Keplerian 
ellipse, by Fourier series. The equations can 
then be treated as a forced harmonic oscillator 
with constant coefficients, and this approximation 
to the equations of motion can be integrated term 
by term. The approximation can be further im- 
proved by using this first-order solution as the 
reference orbit, expanding it in Fourier series, 
and so on. Beyond the second order, however, the 
method is usually replaced by the more general 
one of canonical transformation. 
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The theory of canonical transformation goes 
back to the last century, when it was developed 
to solve mechanics problems . It uses more funda- 
mental properties of dynamical systems, which are 
beyond the scope of this discussion. Suffice it 
to say that it is rigorously equivalent and can 
be either easier or more difficult than the use 
of Lagrange planetary equations. Originally, 
canonical- transformation theory was thought to 
require the use of canonical variables, but the 
generalization by Hori [1966] proved that it can 
be employed for any variable provided the basic 
equations can be solved. In the theory of canon- 
ical transformation, the key is to find a solution 
to a single partial-differential equation. Al- 
though in general this is difficult, a method has 
been developed for the satellite main problem 
that will automatically find a suitable approxima- 
tion to the equation. Therefore, this method is 
applicable to obtaining a solution to any order 
and can be automated, to some extent, on a 
computer. Such an approach is in general use 
now for higher order solutions . 

When the first two methods are inadequate for 
some reason, a third one must be used. This is 
called a semianalytical solution, in that part of 
the solution (integration) can be accomplished 
analytically, while the remaining part must 
necessarily be done numerically. Recourse to this 
method is required when closed-form expressions 
for the force function are complicated or im- 
possible to find. Examples are lunisolar pertur- 
bations when the analytical description of the 
moon's motion to suitable accuracy would be pro- 
hibitive and radiation-pressure perturbations 
when the shadow function must be obtained 
numerically. This method integrates over the 
short-period perturbation analytically and then 
integrates the averaged force function numerically 
to obtain the perturbations. In this case, the 
numerical integration can take relatively large 
time steps and is therefore economical. It also 
conveniently separates long-period and short- 
period perturbations. 

Current Status 

A third-order solution to the main problem — 
i.e., for the motion of a satellite in the 
geopotential containing only Jj, J 3 , and J 4 — has 
been obtained by Kinoshita [1977]". Third-order 
periodic perturbations with fourth-order secular 
perturbations are derived by the method due to 
Hori [1966] . All quantities are expanded into 
power series in the eccentricity, but the solution 
is closed with respect to inclination. A compari- 
son with results obtained by numerical integration 
of the equations of motion indicates that the 
solution can predict the position of a close-earth 
satellite with an accuracy of better than 1 cm 
over a period of 1 month. For this check, a 
special-purpose Taylor-type integrator is adopted, 
in which the positions and velocities are expanded 
into a power series of time and the coefficients 
of the series are determined by recurrence 
formulas [Rabe, 1961; Deprit and Zahar, 1966], 

Periodic perturbations due to tesseral harmonics 
are a first-order linear theory based on integra- 
tion of the Lagrange planetary equations as 


developed by Kaula [1966] . The theory also in- 
cludes the interaction with J 2 and second-order 
interactions with the mean motion through Kepler's 
third law. Although the theory is essentially 
that of Kaula, the details of the calculation 
have been considerably revised with the inclina- 
tion function as described by Gaposchkin [1973] 
or Kinoshita [1977] and the eccentricity function 
calculated as Hansen coefficients. 

The lunisolar perturbations in satellite motion 
are obtained by a semianalytical method [Kozai, 
1973] . The disturbing function is expressed by 
the orbital elements of the satellite and the geo- 
centric polar coordinates of the moon and the sun. 
These coordinates are obtained by using the larger 
terms in Brown's theory [United States Naval 
Observatory, 1954] : 26 terms in longitude, 14 in 
latitude, and 12 in the parallax. The secular 
and long-period perturbations are derived by 
numerical integration, and the short-period per- 
turbations, analytically. Perturbations due to the 
solid body tide can be included in the same way. 

The orbital elements of a close-earth satellite 
have perturbations caused by the motion of the 
equatorial plane of the earth due to precession 
and nutation. Kozai and Kinoshita [1973] derived 
exact differential equations for the perturbations 
of satellite orbital elements due to the motion 
of the earth's equatorial plane and solved them 
to second order in precession. This theory, in 
fact, defines the reference system used for satel- 
lite motion, in which the inclination and the 
argument of perigee are referred to the equator 
of date and the longitude of the ascending node 
is measured from a fixed point along a fixed 
plane and then along the equator of date. 

The perturbations of a spherical satellite due 
to direct solar radiation are computed according 
to a semianalytical algorithm due to Aksnes [1976] , 
which is based on expressions derived by Kozai 
[1961], Through some simple modifications, the 
algorithm also holds when e = 0 and i = 0. The 
perturbations are obtained by summing over the 
sunlit segment of the satellite's orbit during 
each revolution or partial revolution. The end 
points of the segment are evaluated numerically 
once per revolution. Testing of the algorithm 
is done by means of numerical integration of the 
equations of motion and through comparisons with 
observations of the balloon satellite 1963 30D 
during a 200 -day interval. 

The perturbations due to solar radiation dif- 
fusely reflected from the earth have been treated 
by Lautman [1977a] , who used a semianalytical 
method based on the assumptions that the satel- 
lite is spherically symmetric and that solar 
radiation is reflected from the earth according to 
Lambert's law with uniform albedo. Expressions for 
the radiation-pressure force are developed into 
series in true anomaly v. The perturbations with- 
in a given revolution are obtained analytically by 
integrating with respect to v, while holding all 
slowly varying quantities constant. The long- 
term perturbations are then obtained by summing 
the net perturbations at the end of each revolu- 
tion. This theory has been extended [Lautman, 
1977b] to account for the increasing reflectivity 
of the earth toward the poles; the earth's albedo 
is assumed to have a latitude dependence given by 
a = a Q + a 2 sin 2 <p. 
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The short-period perturbations due to air drag 
are computed according to a computation due to 
Sehnal and Mills [1966] . The density function of 
the earth's atmosphere includes the effect of the 
atmospheric bulge described by Jacchia. The . 
method of solution is, in essence, numerical 
construction of the disturbing function. The 
secular acceleration for geodetic satellites is 
more accurately given by analysis of the data. 

The short-period perturbation is usually less them 
1 cm per revolution, and this development is not 
currently used. 

The analysis of ocean tides is done along the 
following lines. Recall that in calculating luni- 
solar perturbations, the body- tide potential is 
easily included by introducing the Love number . 
An alternative way to describe the lunisolar gravi- 
tational potential is essentially given by 
Doodson [1921] . In this case, the potential at 
the earth's surface can be written 


U 

9 


mEp, d, (l)*« ■ 

. Jim tms 
Urns 

__ i(a (t)+mA) 

x P to (*> e 3 


2 

where Re {x} is the real part of x, i = -1, 

[x] is the integer part of x, F^ are numerical 
factors, Djnjs are the s coefficients as determined 
by Doodson by Fourier analysis , P^ m are fully 
normalized associated Legendre functions, and 
a s (t) is the time variation given in terms of the 
six chosen variables. 


Cartwright and Edden [1973] have provided 
values of Dj^s based on modern values of the solar 
and lunar ephemerides: F 20 = -12.020364 cm, 

F 2 i = F^2 = 13.879920 cm. Expression (1) is the 
potential at the surface of the earth, and we can 
continue this potential analytically to satellite 
altitudes as 


U 

g 



^i+m ^ [U+m+l)/2] 


_ i(o (t)+m\) 
x P to( *) e 3 

and then develop satellite perturbations due to 
the sun and the moon. The body tide can be de- 
fined in terms of the complex Love number 



[Munk and MacDonald, 1960, p. 153] . Considering 
the deformation, U t ^ de can be analytically con- 
tinued to satellite altitudes as 


V / a e\^ +1 * 

U tide = g Re f (TV \ F im %s 
Jims ' ' 

x ( i) l+m (-1) I U+m+l)/2] - 

Jim 

i(a ( t) +mA ) 
x e s 

Similarly, each component of an ocean tide of 
height f; s (for driving function D^ ms for argument 
a s ) can be expressed in spherical harmonics ^ ? j r!as / 


where i s fully normalized complex 

representation of the ocean tide. If the ocean 
tide is expressed as a surface layer, then the 
external potential at satellite heights can be 
written 

V 1 +k W a eY +1 

U = 4ttGp a Re L I— J 

ocean me 2Jt + 1 \ r / 


x <i)* +m ( _ 1} IC*-Hn*l)/2J 


x P 4m (40 e 


i(a (t)+mA) 
s 


Jims 


( 1 ) 


As noted by many [Gaposchkin, 1973; Lambeck 
et al . , 1974; Felsentreger et al . , 1976; Goad and 
Douglas, 1978], the ocean and body tides enter the 
potential in exactly the same form, and by satel- 
lite analysis, we can sense only the linear combi- 
nations , 

(s) i + *; (s) _ 

gF im °Jlms k Jl + 4ltGp oj a e 2JI + 1 Jims , 

where k^ and 7 are assumed to depend on the 
frequency of. the argument o s . Thus, we are 
obliged to assume one tide to determine the other. 

For frequencies far from resonance (23* 1 53 m 0?04) , 
k£ and presumably k^ can be taken from seismic 
models [Longman, 1962, 1963; Farrell, 1972] . 

Both Jeffreys and Vicente [1957a, b] and Molodensky 
[1961] pointed out that nearly diurnal earth tides 
should be amplified because of the existence of a 
resonance between the elastic mantle and the 
liquid core. In Table 1, the variation in k 2 
predicted by Molodensky's Model II is given. 

In any event, the expressions for the potentials 
[eqs. (2) through (4)] can easily be used to ob- 
tain satellite perturbations, as the arguments 
are given as linear functions of time, and the 
Lagrange planetary equations can be integrated 
directly as a forced harmonic oscillator. 


Accuracy Assessment 


The various perturbation theories described 
above are all included in a general-purpose 
orbit- determination program that accepts observa- 
tions of direction, range difference, and range 
rate. Each individual theory has been tested, 
and now we wish to study the accuracy of the 
combined theory. For this purpose, we used a 
numerical integration program to calculate simu- 
lated (errorless) data. The general-purpose 
integration program, developed by Krogh [1973], 
is an Adams- type integrator that has the option 
of variable or fixed step size once the integra- 
tion has been started. The variable-step-size 
option uses a desired accuracy as input. The 
program performs the integration in coordinates 
(x, y, z, x, y, z) in single precision on a CDC 
6000 computer that has 14-decimal-digit accuracy. 
The force package allows use of an arbitrary 
gravity field represented in spherical harmonics, 
moon and sun positions (we use the same routines 
as the analytical theory) , radiation pressure, 
and a drag model based on the same physical as- 
sumptions as the analytical developments described 
earlier. 
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TABLE 1. Values of k 2 for Tidal Terms Near Resonance 
(Based on Molodensky, 1961) 


Doodson 

Darwin 

D 

fms 

6 

k 2 

0 (°/hr) 

255.555 

M2 

0.90809 

0.84 

0.3015 

28.984104 

185.555 

00 . 

-0.01624 

-4.54 

0.3025 

16.139102 

167.555 

+1 

-0.00755 

-124.06 

0.3277 

15.123206 

166.554 

n 

-0.00422 

-728.10 

0.4551 

15.082135 

165.565 


-0.0718b 

206.19 

0.2580 

15.043275 

165.555 

K 1 

-0.53011 

192.30 

0 . 2609 

15.041069 

165.545 


0.01051 

180.17 

0.2635 

15.038862 

163.555 

PI 

0.17543 

55.49 

0.2898 

14.958931 

145.555 

01 

0.37694 

6.79 

0.3001 

13.943036 

056.554 

SA 

0.01156 

1.70 

0.3015 

0.041067 



41.87 





y 0. 

2136 - 100 [(a 

+ w)/a] + 

1.9 , 



k 2 = 0.3015 + 0.2109 x 10 3 * 2w) g 

w 


where a = 0 - w, a being the external driving frequency and 6 
being the earth-fixed driving frequency. 


We originally planned to use simulated data 
to verify the short-period perturbations, largely 
because we believed that the accuracy of a 
general-purpose numerical integration could not be 
relied on, but secondarily to conserve computer 
time. For the short-period perturbations, 
simulated range observations were computed for 
subsets of the forces. The simulated observations 
were used in the general orbit-determination 
program, and the mean elements were computed by 
least squares, thus avoiding the problem of 
explicitly relating the initial conditions of the 
numerical integration and the mean elements. This 
testing was done for satellites in orbits similar 
to those of Geos 1 and Geos 3. Our main interest 
was the difference in eccentricity between' the two 
satellites, as one of our concerns was the 
validity of solutions for small eccentricity. The 
results of this testing are given in Table 2 . 

During the first phase of the study, the vari- 
able step size was exercised in the numerical 
integration; the accuracy sought was 10 
cm/sec. Because of the excellent agreement 
of the theory for the main problem (J 2 , J 3 , J 4 ) , 
the integration ephemeris was thought to be of 
sufficient accuracy. However, when working with 
a combined tesseral- and zonal-harmonics field, 
the accuracy never was better than 19 cm root 
mean square. The tesseral-harmonics theory alone 
(with J 2 = 0 ) gave perfect agreement. 

Since their amplitudes are approximately 1 m, 
the J 2 interaction terms could be suspect, but 
they have been carefully checked by me and by 
H. Kinoshita. Furthermore, if the interaction 
terms need revision, we would expect that adding 
more harmonics would increase the error. In fact, 
increasing the field_from_C 2 2 » S 22 to the com- 
plete field through C 44 , S 44 (i.e., eight times 
the number of coefficients) only decreased the 
goodness of fit, from 19 to 22 cm. 

At that point, the accuracy of the integrator 


was questioned, and some simple tests were made 
with fixed-step-size runs. In all cases, the 
trajectories differed by more than 1 m. For 
example, the difference between the variable 
step size and the fixed 0.25-sec step size can 
be written rather well as 

TABLE 2. Test of Orbit Determination 


Force Function 

Interval 

(days) 

Geos 1 
a (m) 

Geos 3 
CT (m) 

J 2 

1 

0.01 

0.01 

J 2 , J3, J4 

1 

0.01 

0.01 

Jv J4 / J5 

1 

0.04 


J 2 ' J 3 ' J 4 

12 


0.06 

J 2' c 22' S 22 

1 

0.19 


c 22' ^22 

1 

0.01 


^ 2 ' ^ 3 * 

C 22 _ ^44 

1 

0 . 22 * 


J 2 » ^3t J 4 * 
^22 - s 44 

6 

0.36 


J 2 + sun + moon 
+ body tide 

1 + 

0.12 


Sun + moon + body 
tide 

1 + 

0.12 


J 2 - ^44 

1 


J 2 ~ ^44 

1 

0.24 



* The integrator used a variable step size. 

T The dynamical effects of the moving equator are 
not included in either the numerical integrator 
or the analytical theory. 

I The integrator used a fixed step size of 0.5 
sec. 

** The integrator used a fixed step size of 0.25 
sec. 


20.1 



e.= 250 an . , 

where t is in days. Two trajectories that differ 
by more than 2.5 m can be fit to an accuracy of 
22 cm. Clearly, the process of fitting mean 
elements was able to adjust or compensate for 
the difference between the two trajectories. This 
difference is due to errors in the numerical- 
integration algorithm, and at this point, it is 
not known what is the true trajectory. The nature 
of the integration error is that much of it can 
be absorbed in the mean elements of the analytical 
theory. In fact, this should not be surprising. 

An along-track acceleration can be modeled by a 
change in the mean motion. In any event, what 
can be said is that short-period perturbations 
can be modeled with an accuracy of 22 cm. The 
limiting factor in this assessment may be either 
the theoretical expressions or the trajectory; 
further investigation is warranted. 

A similar situation obtains in the case of 
lunisolar perturbations, where the principal 
effects are of much longer period than 1 day and 
cannot be evaluated with such a short interval. 

The short-period terms have an amplitude of about 
1 m for Geos-type satellites . From the tests out- 
lined in Table 2, these terms are computed to 
12 cm, or about 10%. We could expect to have more 
than two digits even from the simplest theory. 

The lunisolar perturbations and the tesseral 
harmonics share the common factor that the dis- 
turbing function explicity contains the time. 

This is not true for the zonal harmonics, where 
the disturbing function depends only on position. 
This factor may limit the accuracy of either the 
numerical integration or the analytical solutions. 

Analysis of the Data 

Analysis of tracking data is performed prima- 
rily to obtain geophysical information; a 
secondary consideration is verification of the 
models. (In the latter, I do not include the 
determination of numerical parameters, which fall 
under the primary objective. ) In the results 
given here, both considerations are important. 

The first set of data will concern the Lageos 
satellite, which has exhibited some small and un- 
expected orbital changes . The second set of data 
is on the Geos 1 and Geos 2 satellites, which 
provide useful information about ocean tides and 
core-mantle resonance. Since in both cases the 
results cannot be unambiguously checked, they are 
open to interpretation and are potentially sub- 
ject to errors or oversights in the very complex 
software packages used in analyzing tracking data. 

The Lageos satellite was designed to be a 
stable platform, with a well-defined orbit. Its 
orbital characteristics, given in Table 3, were 
carefully chosen in several respects: The orbit 
is sufficiently high that the effects of the 
anomalous gravity field are reduced, minimizing 
the uncertainty in ephemeris calculation. The 
mean motion minimizes any resonances with the 
gravity field. The very small area-to-mass ratio 
reduces the size of nongravitational perturbations 
(solar pressure, albedo pressure, and atmospheric 
drag) . The symmetrical cross section allows 


TABLE 3. Orbital Characteristics of Lageos 


a = 12.270 x 10 6 m 

e = 0.0046 

I = 109? 86 

n = 6.3866 rev/day 

A/m = 0.0069 cm 2 /g 


greatest ease in attempting to model these forces. 
Most important is that Lageos is equipped with 
cube-corner reflectors, which enables precision 
laser ranging to be done. Orbits for the first 
586 days of the satellite’s lifetime have been 
computed, from which it can be immediately seen 
that indeed the overall objectives have been met 
and the orbit is known very well. 

Figure 1 is a plot of the semimajor axis of 
Lageos . These are independently determined values 
each based on 8 days of tracking data. The per- 
turbations described above have all been included, 
and the remarkable 50-cm decrease in the semi- 
major axis is an unmodeled effect. A candidate 
source for this effect is the use of an in- 
appropriate value of A/m in calculating radiation- 
pressure and albedo perturbations. Over this 
interval, the radiation pressure contributed a 
20-cm decrease in the semimajor axis, and the 
albedo, an additional 10 cm. To compensate for 
the full 50-cm decrease, however. A/m would have 
to be increased by a factor of 2.7, far outside 
the plausable uncertainty in A/m; in addition, 
increases in the specular diffuse coefficient, 
which ranges from 1 to 1.44, and in the solar 
constant would be necessary. A second candidate 
is drag from the. neutral atmosphere. The' 
equivalent drag is enormous, amounting to that 
occurring at altitudes of 2000 km, but adding 
significant drag will reduce the good agreement 
in other orbital elements. In fact, most atmos- 
pheric models do not attempt to model drag above 
2000 km, and our knowledge of atmospheric drag at 
6000- km altitudes is extremely limited. Some sort 
of charge buildup and interaction with the 
magnetosphere is also possible, and a model of 
that should be explored. However, we would expect 
to see some change in a corresponding to changing 
magnetospheric conditions. Nonisotropic radiation 
of heat is another possibility. Again, we would 
expect to see a change as the spacecraft spins 
down. A final possibility is the Poynting 
Robertson effect [Robertson, 1937] . This effect, 
which is viewed in celestial mechanics as an 
aberration, is due to the conservation of momentum 
of photons reradiated from the satellite both 
along and opposite the motion. The effect is 
always along track and has the same effect as 
drag; it can be calculated based on the incoming 
flux or in terms of the temperature at which the 
photons are reradiated. This distinction is 
important if radiation' from the earth contributes 
significantly to the temperature of the satellite. 

Figure 2 presents the variation in inclination 
for Lageos over the same period. Owing to the 
essentially equatorial distribution of laser 
tracking stations, the satellite is not observed 
at maximum latitude and the inclination is not 
so well determined as possible. With improved 
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Fig. 1. Semimajor axis of Lageos. 
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Fig. 2. Inclination of Lageos. 


models, some coming from this work, the inclina- 
tion and other orbital elements will be computed 
with greater accuracy. However, even now, some 
useful information can be derived on the tides; 
this is discussed in the following. 

The study of tides is one of the fruitful 
applications of long-period perturbations . As 


Lambeck et al . [1974] showed, ocean tides as well 
as body tides give rise to sensible perturbations 
in close-earth satellites. Currently, not enough 
information exists about the tides to make a 
definitive analysis, and preliminary results are 
given here to illustrate some of the difficulties 
Though sensible, the perturbations are small and 




TABLE 4. Data Analyzed 



Geos 1 

Geos 2 

Inclination 

59? 383 

105?804 

w 

0?652/day 

-.l?620/day 

u 

-2?246/day 

l?400/day 

First time 

41501 (MJD) 

41501 (MJD) 

Last time 

42239 (MJD) 

42278 (MJD) 

Interval 

738 days 

777 days 


a long time series is needed even to hope to get 
good numerical values. Other physical effects, 
at or near the tidal perturbation frequencies, 
also come into play. For example, some effects 
basically cannot be separated; also, there are 
some very fundamental gaps in our understanding . 

Two series of orbits have been computed on the 
Geos 1 and Geos 2 satellites, as summarized in 
Table 4. in analyzing the mean elements, a 
gravity field and a Love number k 2 are adopted, 
among other constants. We took k -2 = 0.29. The 
inclination of these series has been analyzed and 
interpreted as a perturbation due to ocean tides . 

Individual tidal constituents can be isolated 
by the frequency. In performing harmonic analysis, 
we should also determine terms that should exist 
owing to other effects. For example, odd zonal 
harmonics of the gravity field will not be known 
without error, and hence a term in sin oi will 
arise in the data. Such a term could obviously 
be used to improve knowledge of the odd zonal 
harmonics. Also, we suspect that any analytical 
theory will have an error term that is called 
mixed secular, i.e., t cos ait, and such a term 
should be introduced. 

The harmonic analysis was performed for the 
tidal constituents K2S + K2M, S2, K1S + KIM, and 
PI. The M2 and 01 tides could also be studied, 
but their periods are much shorter, 10 to 15 days. 
Table 5 gives the results of applying equation 
(1) to the observed variations in inclination. 

Also listed in Table 5 are the periods of the 


satellite perturbations, values for some constit- 
uents from terrestrial observations (Lambeck, 
1975], and results of similar analysis by 
Felsentreger et al . [1976] . The fully normalized 
coefficients in complex form are given, together 
with the equi valent amplitude and phase, as is 
conventionally published for ease of comparison. 

In general, the complex tidal component is 
determined to about 1%. 

The first thing to note is that for Geos 2, the 
rates of perigee cl and node i'l are almost equal and 
opposite. Therefore, with a relatively short 
time series , the terms in sin/cos R arising from 
K1S + KIM cannot easily be separated from the 
sin a) term arising from the uncertainty in the odd 
zonal harmonic. The first values in Table 5 for 
that constituent represent an attempt to obtain 
both quantities. The second line corresponds to 
not including a sin w term in the harmonic 
analysis; this result is more plausible and agrees 
with Felsentreger et al . 

A similar circumstance occurs with K2S + K2M, 
in which case, we adopt values not including a 
cos 2oj term, which is equivalent to assuming no 
error in the even zonal harmonics. 

The next point concerns the adopted value for 
the body tide, expressed in terms of the Love 
number Table 6 provides the change in the 

ocean tide's complex amplitude (real part only) 
corresponding to the change in Love number k 2 
from 0.29 (as adopted) to 0.302, the value 
obtained from Farrell [1972], These changes are 
comparable to the derived amplitudes and would 
have a significant effect on any interpretation 
made of the ocean-tide values. Furthermore, if 
we adopt Molodensky’s theoretical model for the 
change of Love number k2 with frequency, then the 
amplitudes of the derived ocean tide would be 
modified by the values given in the last column 
of Table 6. 

Another aspect of the ocean tide is that each 
tidal constituent a s gives rise to a complete 
spherical-harmonics description of the ocean tide 
and a set of harmonics • The satellite orbit 


TABLE 5. Ocean Tides Determined from Satellite Data 


Tide 



Geos 1 




Geos 2 




Lageos 



Darwin Doodson 

Period 

(days) 

(cm) 

c 

(an) 

+ 

E tm 

(deg) 

Period 

(days) 

(cm) 

(cm) 

+ 

C to 

(deg) 

Period 

(days) 

tm 

(cm) 

+ 

(cm) 

+ 

C lm 

(deg) 

K2S + K2H 

275555 

60 

1.711 + 0.931i 

. 1.26 

29 

126 

0.791 - 2.2351 

1.53 

289 












0.598 + 0.5361 

0.52 

42 





S2 

1273555 

56 

2.660 - 0.226i 

1.72 

355 

433 

2.817 + 2.6061 

2.48 

319 

280 

-0.321 + 2.5881 

1.683 

97 

K1S + KIM 

165555 

160 

2.130 + 3.676i 

5.49 

60 

257 

13.260 - 7.3891 














5.30 - 2.441 

7.54 

335 





PI 

1 163555 

85 

-4.779 + 1 . 376i 

6.42 

164 

631 

1.707 +1.1281 

2.63 

33 

221 

-0.349 + 2.644i 

3.43 

98 


Felsentreger 

et al. 

. (1976) 

Lambeck (1975) Geos 1 

Geos 

2 

+ + + + 

C. £. C. e. 

tm tm tm to 

+ 

C ta 

+ 

G tm 




(cm) 

(deg) 

(an) 

(deg) 

(cm) 

(deg) 

S2 

273555 

3.59* 

327 

1.7 

350 

1.0 

62 

K1S + KIM 

1655S5 

3.3 

318 

8.8 

345 

5.7 

334 

PI 

163555 

0.7 

140 

S.O 

178 

4.9 

127 


* Includes atmospheric tide. 
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TABLE 6. Equivalent Ocean Tide 


Tide 

D„ 

4m 

AC 4m (cm) 

k 2 

AC. (cm) 
4m 

M2 

0.90809 

1.908 

0.3015 

1.8288 

K2 

0.11498 

0.242 

0.3015 

0.2316 

S2 

0.42248 

.888 

0.3015 

0.8508 

N1 

-.53011 

-1.114 

.2609 

2.7014 

PI 

.17543 

. 369 

.2898 

-0.0061 

01 

. 37694 

.792 

.3001 

0.6667 


Ak 

= 0.302 - 

0.29 = 0.012 


4 





responds to a linear combination of these ^ , 

where here we have determined the combined effect, 
i.e., a lumped coefficient. To illustrate this, 
consider K2S + K2M. From Geos 1, we have the 
value ~ 1*711 + 0.931i and from Geos 2, 

^p 22 = 0^598 + 0.536i^ If we assume that the only 
other tidal term is ^ 42 * we can obtain the two 
linear equations in two unknowns . It turns out 
that the factor multiplying ^42 for Ge° s 1 is 
negligible, so we immediately have 


<j^? 22 = 1.7ii + 0.931i 


42 = 7 ‘ 44 + 2 * 64i ' 
or in terms of conventional harmonics, 

C 22 = 1*26 cm, e 22 = 28? 56 

C^ 2 = 1.76 cm, e^ 2 = 19?5 . 

These values of C 22 would have to be modified if 
k 2 takes a value different from 0.29. However, 
C^ 2 is unchanged. 

Conclusions 

The accuracy of analytical perturbation theory 
for short-period perturbations seems certain only 
to about 20 cm. It may be better than that, and 
additional, rather straight-forward analysis 
should clarify the situation. In any case, the 
numerical verification of analytical theory is 
not so simple, and the unification of a number 
of different developments based on different 
principals still requires some study. 

The long-term analysis of mean elements is 
providing some new information about the ocean 
and body tides, although many questions remain to 
be resolved. The relation between ocean and body 
tides and the possible frequency dependence of k^ 
and the loaded Love number kjj must be resolved 
before definitive information on ocean tides can 
be obtained. A number of simple tests can be 
made, on which the currently derived values fail. 
They should agree with a numerical model, and the 
relative amplitudes of tides close in frequency 
should correspond to the relative amplitude of 
the driving force as represented by D g . 


Finally, some new information should come from 
analysis of Lageos data. The unexplained secular 
decrease in semimajor axis will not hinder the 
planned purposes of Lageos. In all other re- 
spects, the satellite is being used as it was 
intended. 
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